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L INTRODUCTION 


A 


This  work  was  undertaken  with  the  objective  of  benchmarking  various  turbulence  codes 
against  each  other,  and  to  study  their  sensitivity  to  resolution  and  numerical  techniques. 

The  benchmarking  process  usually  involves  comparison  of  numerical  simulation  results 
obtained  for  the  same  physical  experiment  but  employing  different  numerical  methods  and 
turbulence  closure  approximations. 

In  this  study  the  benchmarking  and  sensitivity  tests  were  performed  on  a test  case 
devised  by  Lewellen  (1974),  which  consists  of  the  collapse  of  a turbulent  mixed  region  In  a 
stratified  fluid.  We  shall  refer  to  this  case  hereafter  as  the  LT  case.  In  particular,  the 
results  of  four  numerical  experiments  will  be  presented:  (a)  a study  of  the  effect  of  mesh 
size  and  stretching  on  the  collapse  process  and  Internal  wave  propagation;  (b)  depth  depend- 
ence resulting  from  the  presence  of  the  free  top  boundary;  (c)  comparison  of  Launder' s 
(1975)  second  order  closure  scheme  with  its  contracted  first  order  version;  and  (d)  com- 
parison of  Lewellen's  numerics  with  ours  on  his  test  problem. 

The  aim  of  experiment  (a)  was  to  find  the  minimum  resolution  necessary  to  simulate 
correctly  the  hydrodyuamic  processes  near  the  collapse  region,  and  the  internal  waves 
propagating  between  this  region  and  the  surface.  This  helped  to  determine  the  smallest 
amount  of  computer  time  and  storage  that  one  could  "get  away  with"  in  running  production 
jobs. 

Experiment  (b)  was  designed  to  study  the  effect  on  the  maximum  velocities  of  the 
presence  of  a free  surface,  serving  as  a reflector  of  internal  waves. 

The  aim  of  experiment  (c)  was  to  establish  the  difference  in  the  generated  internal  wave 
field  when  collapse  is  modeled  by  first-  and  second-order  closure  turbulence  models, 
respectively.  It  also  tried  to  ascertain  whether  first-order  models,  which  require  much 
less  computer  time  and  storage,  would  in  fact  suffice  for  the  problems  of  Interest. 

Experiment  (d),  tried  to  determine  whether  differences  resulting  from  the  employment 
of  various  numerical  schemes  would  in  fact  be  comparable  to  differences  resulting  from 
turbulence  models,  resolution  changes  or  depth  dependence. 

D.  MODEL  DESCRIPTION 

The  lnviscid  fluid  code  used  in  experiment  (b)  Is  identical  to  the  mean  field  part  of  the 

turbulence  code  used  in  the  other  experiments,  with  the  Reynolds's  stresses  set  equal  to 

zero.  These  equations,  together  with  the  corresponding  numerical  schemes,  have  been 
Note:  RSnuscript  submitted  December  11,  1978. 


presented  In  a paper  by  Dugan  and  Wam-Vamas  (1974)  and  will  not  be  discussed  separately 
here. 

Our  second-order  closure  approximations  for  turbulence  modeling  follow  the  Launder 
(1974)  philosophy.  The  main  difference  between  this  scheme  and  that  presented  by  Lewellen 
(1974)  is  in  the  closing  of  the  triple  velocity  correlations,  and  the  correlations  between  the 
pressure  and  the  velocity  gradients  or  the  velocities. 

We  shall  now  proceed  to  describe  the  Launder  turbulence  model.  We  partition  all 
dependent  variables  Into  their  mean  and  fluctuating  parts,  e.  g.  m = U(  + U[*,  and  write 
transport  equations  for  the  time  evolution  of  the  mean  fields  and  the  second-order  correla- 
tions, e.  g.  uiuj : 

Dffi  1 3p  P 3 

+ — gj (ujuj)  + vV^uj  l1 

Dt  P0  3xt  p0  3Xj 


N Triple  velocity  correlation! 


ulwj“k  " " cl7^»iuk  ^ <uiuj).  «|  • .U 


set  equal  to  zero.  The  neglect  of  pressure  diffusion  has  some  experimental  Justification, 
with  the  most  persuasive  being  that  of  Irwin  (1974). 

The  length  scale  of  the  energy  containing  eddies  is  obtained  from  a relation  between  the 
turbulence  energy  density  K and  the  dissipation  function  c : 


t ■ 


rS/2 


(14) 


Lewellen  (1974)  makes  a somewhat  different  closure  approximation  for.the  correlations 
f Cl  1 - I C4  1,  using  the  relation  (14).  For  the  pressure-strain  correlation  [ Cl  ] he  sets 
c2  = 0 in  (6),  which  essentially  omits  terms  that  represent  the  effects  of  gravity  and 
production  of  turbulence  by  the  mean  velocity  shear.  A similar  effect  is  achieved  by  setting 
c2t  = 0 in  (10)  for  the  pressure-density  gradient  correlation  [ C3  ).  In  the  correlations 
[ C2  ) and  [ C4  ) Lewellen  neglects  the  terms  involving  u£u£  . where  1 t k,  and  for  the 
remaining  term  £ u£u£  he  substitutes  from  (7a)  and  (14).  Thus,  ( C2  | becomes 


ulujuk  " ~ 2 c3  *K  '2  3^  <uiuP 


(15) 


and  ( C4  ) 


■ -2c 


3t 


£K 
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(16) 


To  obtain  a first-order  closure  scheme,  we  will  make  the  following  closures  of  the  second- 
order  correlations: 


— — v2  /^Ui  \ 

uluj  • - cu  — y3Xj"  + 3 ^ j 


(17) 


Cu  K2  dp 


u.p 


— - - — « c„  ■ .09  and  o.  • 1.0 

c 3xj  * u 0 


(18) 


The  predictions  of  u '2  proceeds  in  the  same  way  as  in  the  second-order  theory. 
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ID.  RESULTS  AND  DISCUSSION 


A.  Mesh-size  Dependence  Maximum  Velocities 

The  second-order  closure  model  of  Launder  (1975)  described  above  was  used  to  study 
the  sensitivity  to  resolution  of  the  wake  region.  In  order  to  save  computer  time.  Initial 

r 

conditions  somewhat  different  from  Lewellen  (1974)  were  employed. 

The  computational  region  consisted  of  a square  of  200  cm  by  200  cm.  with  the  radius  of 
the  mixed  region  being  a = 1 5 cm.  A linear  stratification  of  N - . 707  sec  was  specified 
everywhere  In  the  fluid  at  t = 0,  and  an  Initial  patch  of  turbulence  energy  was  deposited 
according  to  the  distribution  law 


(19) 


where  A ” 103cm-/sec2  and  a is  the  radius.  The  mean  velocities  and  til  second-order 

correlations  were  set  to  zero  at  t = 0. 

Figures  2 and  3 show  the  comparison  of  the  maximum  horizontal  and  vertical  velocities 

obtained  with  resolutions  of  5,  7.  and  9 grid  points  in  the  (Initial)  wake  radius.  We 

observe  an  approach  to  asymptotic  stabilization  of  the  results  for  7 or  more  mesh  points. 

The  percentage  deviation  of  the  5 grid  point  results  from  the  asymptotic  results  is  »%  for 

both  u and  w 

max  max 


B.  Depth  Dependence  of  Maximum  Velocities 

This  experiment  was  carried  out  by  the  invlscid  code  reported  on  by  Dugan  and  Warn- 
Vamas  (1974)  In  NRL  Memo  Report  #2597.  The  resolution  of  the  wake  was  8 grid  points  ir. 
the  vertical  and  7 in  the  horizontal  in  all  cases.  Figures  4 and  5 show  the  results  for  depth- 
to  radius  ratios  of  2 and  8,  respectively;  the  percentage  variation  in  the  maximum  velocities 
is  about  8%.  These  values  occurred  near  the  edge  of  the  wake,  and  can  be  associated  with 
the  circulation  that  takes  place  during  the  primary  collapse  phase  of  the  wake. 

It  was  also  Important  to  determine  whether  the  results  would  change  further  if  the  depth- 
to-radlus  ratio  of  8 is  increased.  To  simulate  the  asymptotic  case  of  an  infinite  ratio,  an 
absorbing  liner  was  placed  outside  the  wake;  mathematically  this  is  represented  by  a region 


of  high  Newtonian  viscosity.  The  following  expression,  suggested  by  Roberts  (1972),  was 

" - *5  exp j SINH  ^SINH  (-  x/.125  Lx)j  + S1NH  J^SINH  (-s/.125  Lf)J  j (20) 

where  L and  L are  the  dimensions  of  the  computational  domain  in  the  x and  z regions, 

X E 

respectively.  The  results  were  found  to  be  equivalent  to  those  of  the  ratio  of  8. 


st  nd 

C.  Comparison  of  1 and  2 - Order  Closure  Results 

These  experiments  were  performed  with  the  launder  (1975)  turbulence  model  described 
in  Section  II  above.  The  test  case  on  which  they  were  compared  was  devised  by  Lewellen 
(1974),  and  is  referred  to  in  this  paper  as  the  LT  case.  The  computational  domain  was  a 
square  region  of  area  10“  X 10“  cm-.  Because  of  the  inherent  symmetry  of  the  problem  in 
both  x and  z in  the  case  of  a linear  ambient  density  stratification,  only  a quarter  of  the 
whole  region  needs  to  be  simulated.  The  "radius"  a of  the  turbulent  region,  i.e.  the  half- 
width of  a gaussian  distribution  of  perturbation  density,  was  500  cm,  making  the  depth-to- 
radius  ratio  of  this  problem  20.  A scaling  analysis  of  Eqs.  (1)  - (18)  was  performed  using 


Here  u is  the  free 


a length  scale  a,  a velocity  scale  u and  a density  variation  scale  a ' o 

Je~ 

stream  velocity  of  a self-propelled  object  that  in  fact  never  appears  explicitly  in  the  model, 

using  a time  scale  t ■ a/u  the  following  nondimensional  parameters  emerge: 
o 


N2 


8 3P° 


FT  - 


constant. 


(21a) 

(21b) 


aN 


Ri 


Fr**»ax 


(21c) 


Here  K'  is  the  maximum  value  of  the  initial  nondimensionalized  turbulence  energy  den- 
max 

sity,  N is  the  Brunt- Vaisaln  frequency,  Fr  is  the  Froude  number  and  Ri  is  a Richardson 

number  based  on  both  the  stream  velocity  and  the  turbulence  velocity.  In  nondimensional 

form,  the  initial  conditions  become 

5 ■ z * 5.  1 


K - f (.0108)  e 
A 


-2(r2  - 1) 

- 1.38  r2 


r2  - x2  + z2 


(22) 

(23) 
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for  the  perturbation  density  and  the  turbulence  energy,  respectively. 

To  start  a specific  calculation,  we  have  to  specify  Fr  and  Ri,  or  Fr  and  K'max-  10  the 

experiments  performed  Ri  = .872  and  Fr  • 10.5.  All  correlations  of  form  ui'p'  and  u£uj 

with  i * j were  set  zero  at  t = 0.  Finally,  one  has  to  specify  the  integral  length  scale  of 

turbulence  t introduced  in  (14),  and  also  referred  to  as  the  macroscale.  It  is  determined 

by  relating  it  to  a fraction  of  the  distance  from  the  center  of  the  wake  to  where  the  turbulence 

drops  to  a quarter  of  its  maximum  strength.  In  practice,  one  finds  both  the  length  scales 

j and  t that  measure  the  distance  from  the  axis  to  the  isoline  K = .25  K of  turbu- 
x max 

lence  along  the  x and  z directions,  respectively,  and  then  obtains  t according  to  the  formula 


.2  tjl 


xlz 


{x  + i‘z 


(24) 


Figures  7 and  8 show  the  first  and  second  order  closure  results  for  u , w , P 

° max  max  max 

and  o . where  q = (2  K ) S is  the  turbulence  velocity.  The  location  of  the  maxima 
max  max  max 

of  the  velocities  in  time  are  approximately  the  same,  but  the  second  order  magnitudes  are 
somewhat  higher.  The  second  order  density  perturbation  is  also  higher,  whereas  the  turbu- 
lence velocity  is  slightly  lower.  Table  I gives  the  exact  numbers  at  t a .35  nondimensional 
ttme.  Figure  9 illustrates  the  time  development  of  the  various  forms  of  energy  in  time. 

Both  the  kinetic  and  available  potential  energy  associated  with  second  order  closure 
exceed  the  first  order  values  by  almost  a factor  of  2.  The  turbulence  energy  curve  is 
the  same  for  both  closures,  indicating  a possible  dominance  of  the  decay  term  in  the 
turbulence  energy  equation.  The  available  potential  energy  is  defined  as 


where  is  the  Initial  mean  density,  and  the  deviation  from  Px  represents  the  largest 
fraction  of  energy  in  each  model.  A possible  explanation  for  the  potential  energy  behavior 
is  given  further  below,  in  connection  with  an  analysis  of  the  magnitude  of  terms  in  the 
density  Eq.  (2). 

Next  we  proceed  to  a comparison  of  the  l"  * and  2 -order  results  for  the  temperature- 
velocity  correlations  and  the  off-diagonal  Reynold's  stresses.  Figure  10  exhibits  the  turbu- 
lent transport  of  density  fluctuations  v'q'  as  a functlon  z at  */a  = ,57  at  E " 
.37.  Differences  in  both  magnitude  and  shape  are  observable,  with  the  second  order  trans- 
port everywhere  negative  and  the  first-order  reversing  sign  half-way  to  the  vanishing  point. 
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Since  in  1 -ordor  theory  w'p'  - - K ||  , the  turbulent  transport  follows  the  sign  of 
so  that  outside  the  wake,  up  to  r ■{  2 a,  there  is  a counter- gradient  density  flux  according 
to  2 -order  theory.  This  would  imply  that  negative  entrainment  is  being  generated  by  the 
first-order  model,  i.e.  some  kind  of  free-convective  turbulence  that  produces  mean  poten- 
tial energy. 

Figure  11  displays  the  Reynold's  stress at  x/a  = .57,  t'  = 37.  In  this  case  there 

are  noteworthy  differences  in  the  magnitude,  shape  and  sign  of  the  stresses,  with  the 

second-order  stress  being  everywhere  negative  and  the  first-order  positive  out  to  about 

1.3  radius.  However,  since  the  effective  viscous  dissipation  is  given  by  — (u'w')  it  is 

3z 

the  slopes  of  the  two  curves  that  should  be  compared.  They  both  exhibit  a negative  slope 
out  to  about  z ^ a,  i.e.  the  edge  of  the  (original)  wake,  indicating  normal  viscous  dissipa- 
tion. Then  for  a ^ z^  1.5a  in  both  cases  the  slopes  become  negative,  Indicating  a pos- 
sible mean  flow  enhancement  by  the  eddies.  The  big  discrepancy  occurs  in  the  magnitude 
of  the  dissipation,  with  the  first-order  results  exceeding  the  second-order  values  by  a 
factor  of  ten.  In  the  case  of  the  buoyancy  flux  ratio  in  Figure  10  this  ratio  is  only  ~ 3. 

Next  we  investigated  the  magnitude  of  the  individual  terms  in  the  mean  density  trans- 
port Eqs.  (1).  The  results  are  summarized  in  Figure  12,  showing  each  term  as  a function 
of  horizontal  distance  at  a height  z/a  = . 587  above  the  center  of  the  wake,  at  time 
f ■ Nt/2TT  » .37.  It  is  evident  that  transport  by  the  'thermal  stresses'  or  'buoyancy 
fluxes  u'p  ' and  w'p 'is  comparable  to  that  by  the  mean  flow  in  the  first-order  results. 

On  the  other  hand,  in  the  second-order  results  the  turbulent  contributions  are  two  orders 
of  magnitude  less  than  the  mean  advection  terms.  This  great  discrepancy  between  first- 
and  second-order  closure  simulation  explains  the  difference  by  a factor  of  two  in  the 
available  potential  energy. 

D.  Evaluation  of  Lewellen  and  Warn-Varnas  Numerics  on  the  LT  Benchmark  Case 

It  was  of  inte  rest  to  know  how  the  present  code  would  perform  on  the  LT  Benchmark 
case  if  we  set  the  turbulence  closure  and  pertinent  constants  identical  to  Lewellen's,  effect- 
ively comparing  the  different  numerics.  To  this  end  we  have  set  the  closure  constants 
c2t  “ c2  = 0 if,  equations  (1)  - (12),  leading  to  an  exclusion  of  gravitational  contributions 
to  the  pressure-strain  correlation.  The  other  constants  have  also  been  made  to  agree  with 
Lewellen's  values.  Using  equation  (14)  relating  the  turbulence  length  scale,  dissipation 
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function  and  turbulence  Intensity,  the  pressure-strain  correlation  (6)  becomes 


L 


where  q = \|2iT 


The  diffusion  term  In  (4)  involving  the  triple  correlation  in  the  reduced  form  given  by 
(15)  simplifies  to  ^39  q*.  j<ti  , with  <t>  being  the  stress  calculated  and  c3  set  to 
.28.  The  dissipation  correlation  for  Lewellen's  case  was  closed  as 


9u£  3uj 

with  a = 2. 5 and  b = . 125. 


2av 

t2 


The  calculation  was  performed  on  a region  20a  by  15a  (a  being  the  wake  radius),  with 
a 40  x 40  grid  resolution  and  placing  8 points  of  a stretched  mesh  inside  the  wake.  The 
results  are  summarized  in  Figure  16  and  Table  I (4th  line).  The  differences  in  the  maxi- 
mum horizontal  and  vertical  velocities  are  8%  and  1 9%,  respectively.  There  is  also  a 
shift  in  the  occurrence  times,  from  t'  = Nt/2*  of  .35  and  .21  to  .37  and  .24,  respectively, 
with  Lewellen's  code  giving  the  later  values. 


Comparing  Figure  16  to  Figure  6,  the  most  striking  difference  is  the  slower  decay  of 
the  turbulence  intensity  q for  the  Lewellen  numerics  (Figure  6).  The  most  logical  source 
of  this  discrepancy  would  be  a difference  in  the  constants  occurring  in  the  diffusion  and 
decay  terms  discussed  above.  The  description  of  the  macroscale  l would  also  influ- 
ence the  decay.  The  differences  in  the  q-decay  would  also  influence  the  behavior  of  the 
maximum  horizontal  and  vertical  velocities. 


rv.  CONCLUSION 

Comparison  calculations  have  been  made  on  the  collapse  motions  and  Internal  waves 
generated  by  turbulent  wakes  deposited  In  a stratified  fluid.  The  main  results  of  the  cal- 
culations can  be  summarized  as  follows 

1.  Seven  mesh  points  in  a wake  radius  are  sufficient  to  resolve  the  collapse  prooess, 
with  five  points  giving  an  8%  error  in  the  maximum  velocities  generated. 

2.  A change  in  the  depth- to- radius  ratio  of  two-to-eight  results  in  an  8%  variation  in 
the  maximum  velocities.  Any  further  increase  In  the  ratio  results  in  no  farther 
veloolty  variations.  These  results  imply  that  the  presence  of  the  top  surfaoe  (a 
density  discontinuity)  has  no  great  Influence  on  the  collapse  process,  provided  the 
maximum  extent  of  a turbulent  wake  is  less  than  half  Its  depth. 

3.  A comparison  of  1st  and  2nd  order  turbulence  closure  schemes  showed  a 33% 
variation  in  the  maximum  horizontal  velocities  and  a 25%  variation  in  the  vertical 
component,  indicating  a need  for  a better  tuned  or  physically  oorrect  lst-order 
model.  Furthermore,  a comparison  of  the  correlations  w~r”  and  u 'v'  with  the 
second-order  results  gives  evidence  of  negative  density  entrainment  (free- convective 
energy  release)  and  mean  flow  enhancement  by  the  Reynolds  stress,  respectively, 
features  that  lst-order  theory  can  not  reproduce.  However,  we  have  not  deter- 
mined the  total  region  size  of  the  Quid  where  such  processes  take  place,  nor  their 
total  quantitative  contribution  to  the  collapse. 

4.  Comparison  of  numerics  by  Wam-Varnas  and  Lewellen  show  a discrepancy  of  20%, 
which  made  the  comparison  of  the  Launder  closure  with  the  Lewellen  turbulenoe 
closure  difficult  to  interpret. 


Table  I 


Comparison  of  maximum  mean  velocities  for  the  LT  case  using  different  closure 
schemes  and  numerics.  The  horizontal  velocity  umgx  «nd  wmgx  are  given  in  di- 
mensionless units,  and  their  respective  deviations  from  the  LT  results  in  %. 


umax 

^max 

«U(X) 

«W(X) 

LT  (Lewellen,  1974) 

.0310 

.0285 

0 

0 

1st  order  closure 

.0260 

.0210 

-16 

-26 

2nd  order  Launder  closure 

.0370 

.0280 

19 

-2 

2nd  order  Lewellen  closure 

.0285 

.0235 

-8 

-18 
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APPENDIX  A 


FINTT E-DIFFERENCE  SCHEME 

We  shall  first  discuss  the  stretching  of  the  mesh  and  the  spatial  differencing,  then  the 

time  differencing  and  the  method  used  to  find  the  pressure.  The  stretching  of  the  mesh  was 

x/x  x/x 

accomplished  by  the  use  of  the  function  £.(e  i-i)/(e  4 + i ) where  x l Is  a 

normalization  oonstant  that  oontrols  the  stretching  In  the  £ -th  boundary  layer.  For  equal 
Increments  At  of  the  function  - iAf.,  a set  of  coordinate  values  xt  will  be  generated  suoh 
that  A x=xj  - x 1_1  constitute  the  variable  mesh  spaolngs.  A similar  transformation  la 
effected  between  n and  z . It  must  be  pointed  out  that  the  equations  are  not  transformed 
to  a new  coordinate  system  £ and  n , but  are  directly  differenced  on  the  new  grids  xj  and 
ZV  ; the  funottons  f.  and  n merely  serve  the  purpose  of  creating  a smoothly  varying 
x and  7 mesh. 

The  differencing  of  the  adveotlon  terms  will  be  Illustrated  on  the  transport  part  of 
the  horizontal  velocity  equation.  The  scheme  was  presented  for  constant  grids  by  Placsek  6 
Williams  (1970).  We  shall  denote  the  coordinate  axes  by  xj  and  zk  , the  time  by  tn=n  A t, 
and  the  value  of  a dependent  variable  $ as  4>  (xj,  zk  5 tn)  = 4,5^  with  the  superscript 
denoting  a time  level  and  the  subscripts  a mesh  point  location  In  space.  Reference  should 
be  made  to  Figure  1 for  inspecting  the  arrangement  of  the  dependent  variable  on  the 
various  subsets  of  the  staggered  mesh.  (Note  that  upper  oase  symbols  are  used  to  denote 
the  variables  In  the  figure. ) Using  these  symbols,  we  have. 


7- — — — [(“i+i  + “i^i+i  - (“i  + “l-l>“i-il  + 

4Ax1+i}  l Jk 


+ + i'lJs)k'H*  “1*k+1  ‘ <i5l44s  + "i.k-l]. 

The  advantage  of  this  scheme  is  that  It  conserves  u exactly. 


E-  ^U1  k 

Ui>  AXlA'k  ’ ° 
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when  summed  over  the  interior  grid  points,  regardless  of  whether  or  not  the  divergence  is 
zero.  For  a discussion  of  this  point  and  the  method  of  evaluating  conservation  properties 
the  reader  is  referred  to  Williams  (1969)  and  Piacsek  and  Williams  (1970).  It  is  known  that 
the  boundedness  of  u 2 ensures  stability  and  that  of  u does  not;  furthermore,  nonzero  diver- 
gence is  always  present  when  iterative  methods  are  used  to  find  the  pressure.  The  remain- 
ing terms  in  the  horizontal  velocity  equation  are  differenced  as  follows : 


the  pressure  term. 


the  Reynolds  stress  term. 


1_ 

^xi+*s 


zk+lj 


[(u'u')i^  - <u^)j:  -I 


the  viscous  term, 


” [i»1+ , <s?<  i.k  ■ “ltk1  - 5^  <“i,k  - “ili.k’J 

* [i.kt , (“?.k..  - - i,k  <“?.k  - ] } 

The  finite-differencing  of  the  forecast  equations  for  the  Reynolds  and  thermal  stresses  is 
illustrated  on  the  sample  equation, 

h + h [*!ru'u'] 

- F u "u ' + R 

in  finite-difference  form  this  equation  becomes: 

u-u'^Vu'*-  1 _ 1 ( 4’1+ls  ”\n  f-r-.n+ll 

-w—  ■i*J^r)i+'‘(uu)i  1 
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+ R 


n 


where  Rn  represents  the  advectlon  and  all  other  terms  evaluated  at  time  level  n. 

The  Poisson  equation  Is  formed  from  the  finite-difference  equations  by  performing  the 
divergence  operation  on  them  In  finite-difference  form  (Williams,  1969).  The  divergence 
is  forced  to  zero  within  round  off  error  by  adjusting  the  pressure  through  the  inclusion  of 
the  divergence  In  the  source  term  at  time  level  n-1. 


The  Poisson  equation  Is  solved  by  an  ADI  iterative  approach  as, 

(r.  - )pJ  k - (r*  + k 


(r  _ - <r0  + w - S.  k 

t u2>  i.k  t 3x2  i.k  i.k 


where  r4  are  the  iteration  parameters  with  t denoting  the  iteration  numbers.  Pi>k  Is  the 
pressure  and  k the  source  term  on  the  grid  point  (i,k).  The  continuum  second  deriva- 
tive operators  are  understood  to  represent  their  appropriate  finite-difference  analogues. 

The  boundary  conditions  on  the  pressure  are  of  the  Neuman  type,-i£  - c,  where  s is  the 

a S 

distance  normal  to  the  boundary  and  G is  obtained  by  applying  the  mean  field  equations  on 
the  boundary.  The  optimum  Iteration  parameters  are  calculated  by  the  method  outlined  in 
Wachpress  (1966).  This  method  minimizes  the  maximum  eigenvalue  of  the  Iteration  matrix. 
The  truncation  error  of  the  overall  finite-difference  scheme  Is  between  first  and  second 
order  in  time  and  space. 

The  restrictions  on  the  time  step  are  derived  empirically  by  trying  different  time  steps 
for  a given  problem  of  a given  turbulent  strength  and  Brunt- Vais sala  frequency.  The  re- 
striction on  the  time  step  for  some  types  of  terms  as  the  non- linear  transport  terms, 
which  are  subject  to  the  Courant-Friederlcks-Lewy  condition  of  AtJ*  A s/vD  where  As  Is 
the  smallest  grid  spacing  and  v0  Is  the  largest  velocity.  Is  known.  Similar  Insight  can  be 
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gained  Into  the  stability  requirements  due  to  some  types  of  coupling  between  the  equations 
that  govern  Internal  gravity  waves,  but  even  a partial  linearized  analysis  of  the  equations 
becomes  rapidly  too  time  consuming.  The  non-physical  computational  mode  that  arises 
from  the  flnlte-dlfferenclng  of  the  non-linear  terms  where  a first  order  equation  in  time 
is  raised  to  a second-order  difference  equation  is  eliminated  by  periodic  averaging  of  the 
variables. 

The  boundary  conditions  on  the  mean  field  velocities  are  that  the  normal  velocity 
vanishes  and  the  normal  derivative  of  the  tangential  velocity  vanishes.  For  the  Reynold 
and  thermal  stresses,  we  use  the  boundary  condition  of  vanishing  derivatives, 
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HORIZONTAL  DISTANCE 


(Nt/27T) 

Fig.  2 - Comparison  of  maximum  mean  horlaontal  and  vertical  velocities 
obtained  with  different  wake  resolutions.  Case  a)  places  a constant 
grid  of  5 points  in  the  wake,  and  case  b)  a stretched  grid  of  9 points. 
The  radius  In  each  case  Is  15  cm  and  the  total  computed  region  Is  100 
cm  x 100  cm. 
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Fig.  3 - Continuation  of  Figure  2,  with  case  a)  again  placing  a 
constant  grid  of  5 points  in  the  wake  and  case  b)  7 points  on  a 
stretched  grid.  Here  case  b)  had  a total  region  sise  of  200  cm 
x 200  cm. 
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MAXIMUM  HORIZONTAL  VELOCITIES  (cm /sec) 


MAXIMUM  VERTICAL  VELOCITIES  (cm /sec) 


VALUES 


ENERGETICS  SCALE 


.1333  .2666  .3999 

(Nt/27r) 


Fig.  9 - Comparison  of  tha  tint  evolution  of  available  potential  energy 
(APE),  kinetic  energy  (KE)  and  turbulence  energy  (TURB)  for  let  and  2nd 
order  turbulence  closure  scheme  a 
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Fig.  10  - Plot  of  the  thermal  strata  w 'p ' va.  depth  at  a horiaontal 
distance  of  x/a  “ .57.  Curve  a rapresanta  lat  order  closure  and  curve 
b 2nd  order  closure. 
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Fig. 


15  - Contours  of  assn  density  (deviation  from  horlsontal  average) 
for  the  LT  case,  obtained  with  second-order  closure 
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Fig.  16  - Plot*  of  maximum  mean  density,  mean  horlsontal  and  vertical 
velocity  and  turbulent  energy  vs.  time  for  the  LT  case,  but  with  the 
authors'  closures  and  constants  set  to  agree  with  Lewellen  (1974) 
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